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In cardiac electrophysiology, it is important to predict the necessary conditions for conduction fail- 
ure, the failure of the cardiac excitation propagation even in the presence of normal excitable tissue, 
in high-dimensional anisotropic space because these conditions may provide feasible mechanisms for 
abnormal excitation propagations such as atrial re-entry and, subsequently, atrial fibrillation even 
without taking into account the time-dependent refractory region. Some conditions of conduction 
failure have been studied for anisotropy or simple curved surfaces, but the general conditions on 
anisotropic curved surfaces (anisotropic and curved surface) remain unknown. To predict and ana- 
lyze conduction failure on anisotropic curved surfaces, a new analytic approach is proposed, called 
the relative acceleration approach borrowed from spacetime physics. Motivated by a discrete model 
of cardiac excitation propagation, this approach is based on the hypothesis that a large relative ac- 
celeration can translate to a dramatic increase in the curvature of the wavefront and, subsequently, 
to conduction failure. For simple anisotropic surfaces, the relative acceleration approach is validated 
by computational simulations or the previously known results from the kinematics approach. As 
a practical application, this approach is proposed to provide theoretical explanations of the mech- 
anism of cardiac excitation propagation around the pulmonary vein with anatomically observed 
anisotropy. 
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I. INTRODUCTION 

In modern cardiology, the structure and myofibers of the atrium are gaining a broader audience. Ho et al. ^ 
studied the different topological structures and the orientation of myofibers of the atria to reduce possible risks 
associated with interventional procedures. Anderson et al. [2 also stressed the importance of the knowledge of the 
atrial structure with anisotropy as it correlates to muscular activity and response to electrical signals. There are 
more specific questions: does anisotropy in the plane play the same role as it does on the surface of the atrium [6] 
[41 ? How is the reentrance, or the reoccurrence of cardiac excitation, generated around the pulmonary veins of the 
unique geometrical shape with specially-oriented myocardial tissue ^S] [34] [39] [44] ? In mathematical physics, these 
problems are all related to the behaviour of cardiac excitation propagation, as a system of diffusion-reaction equations, 
in high-dimensional anisotropic space. 

The mathematical analysis of the wave propagation in excitable media was first obtained by means of the kinematics 
of the wavefront, known as the kinematics approach^ to analyze the effects of the geometrical shape or anisotropy on 
the behaviour of excitation propagation. As first appeared in ref. [12] in English, the kinematic approach has been 
used to study the propagation of excitation in inhomogeneous and anisotropic media. The kinematic equation defines 
the relationship between the curvature of the wavefront and arc length on curved surfaces to yield the conditions of the 
critical curvature which breaks up the propagation of the wavefront on the curved surface. The kinematics approach 
has been successfully used to describe the geometrical effects on surfaces such as non-uniformly curved surfaces [11], 
periodically modulated curved surfaces [13 , ring-shaped propagation on curved surfaces [14 , propagation on a torus 
[15], and on moving excitable media [16]. In addition, this approach was also used to analyze the role of anisotropy 
in the plane: the breakup of the propagation [29 and the propagation of curved fronts in anisotropic excitable media 

However, here some problems begins. For anisotropy on curved surfaces, or when considering both of anisotropy 
and curved surfaces, the kinematics approach becomes too complicated even on a simple surface. Consequently, the 
higher dimensional anisotropic space, which is a more realistic approximation to the atrium and ventricle, seems to 
remain beyond the scope of the kinematics approach. Consider the following kinematic equation [12 to describe the 



2 



relationship between the curvature of the wavefront (K) and arc length {£) on curved surfaces as 

where V is the normal propagation velocity, T is the local Gaussian curvature of the surface, and C is the tangential 
growth velocity of free ends of the front. Even disregarding the anisotropy and the shape of geometry affecting the 
multiple variables in the above equation, the most difficult complexity arises due to the local Gaussian curvature (F) 
which depends on both anisotropy and geometry. This complexity remains the same even for the time-independent 
case such as the spiral wave. To overcome this problem, a new analytical tool has been in demand and this motivates 
this paper. 

A. Prom wavefront to trajectory 

The new approach begins with introducing the concept of the trajectory into cardiac excitation propagation. This 
implies that cardiac excitation propagation will be regarded as a (physical) wave. Briefly stated, the trajectory is 
the path of a particle of the wave which is described as same as in the context of classical mechanics. If we put a 
particle on the wavefront of cardiac action potential at time to, the continuous observation of the particle in time 
t > to generates a trajectory. Consequently, in homogeneous media, the direction of the trajectory is orthogonal to 
the direction of the wavefront. 

The definition of the trajectory will be provided more rigorously in Section 2, but this meaning of the trajectory 
can be ambiguous intrinsically in cardiac excitation propagation. The concept of the trajectory requires the presence 
of traveling particles, but it has been assumed that no particles actually propagate through myocardial tissues on 
macroscopic scale. There are complex movements of several ions, such as A/'a+, Ci~ , Co?^ , on microscopic scale 

through ion channels and gap junctions, but it has not been accepted that certain kinds of ions can be identified as 
travelling particles for cardiac excitation. However, a trajectory can be uniquely constructed by a smooth wavefront, 
which is well defined clinically and mathematically, in a sufficiently smooth surface as shown later in this paper. This 
is another reason we stick to sufficiently smooth surfaces, but the concept of the trajectory still remains the same as 
any other physical wave involving particle motion. 

The relative acceleration approach is therefore based on the study of the trajectory in cardiac excitation propagation. 
The aspect of this approach reminds us that the kinematic approach concerns the shape of the wavefront, but the 
relative acceleration approach concerns the distance between the trajectories. From this perspective, the relative 
acceleration approach can be regarded as a complementary method to the kinematics approach, since the distribution 
of the trajectories reflects the shape of the wavefront. Nevertheless, we will show that the relative acceleration 
approach can be expressed a lot simpler than the kinematic approach, especially for anisotropic curved surfaces and a 
higher-dimensional anisotropic space. This occurs because the trajectory can be easily written in terms of the metric 
tensor which incorporates anisotropy and curved space into the same geometric term. 

B. Regarding a priori properties and restrictions of the surface 

In this paper, the relative acceleration approach will be applied only to anisotropic curved surfaces for the following 
reasons: easy derivation of the governing equation, easy validations of the analysis from the previously known results, 
and computational modelling. The extension to a higher dimensional anisotropic space can easily be derived and 
will be reported in later publications. Note that the choice of the domain as curved surfaces does not necessarily 
mean that the real heart, especially the atrium, can be successfully modelled by curved surfaces. This issue is in 
fact beyond the scope of this paper. Regard this paper as the first step of the proposed approach toward generally 
higher-dimensional anisotropic spaces. 

For the sake of simplicity of analysis, two a priori properties are made on the curved surface M unless mentioned 
otherwise. 

The first a priori property is inferred from an assumption that the derivation of the equations is based on M. which is 
two-dimensional manifold and locally- Euclidean. Roughly stated, the curved surface is a two-dimensional manifold 
in the sense that, for the neighbourhood of each point on the curved surface, there exists a one-to-one representation 
of each curved element onto a simply connected region of the Euclidean plane [7j. Moreover, the curved surface M is 
locally-Euclidean in the sense that the curved element can be represented as a small domain of the Euclidean space 
in a small neighbourhood of any point 0. For example, the surface of revolution is always manifold and locally 
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Euclidean. All these assumptions are satisfied, alternatively, if there is a pair of orthogonal curved axis at every point 
as frequently used in curved spacetime. 

The second a priori property is that all quantities and equations are locally defined and are only piecewisely con- 
tinuous. Since all the derivations and definitions are achieved in a curved element of arbitrary size from a tessellation 
of the curved surface, the corresponding quantities and properties must be well-defined within each element, but do 
not have to be continuous and different iable across the elements. For example, the supposition on the orthogonality 
of the curved axis, the existence of orthogonal diffusivity tensor, and the existence of trajectory for smooth wavefront 
can be easily justified point-wisely due to the fact that the surface is assumed to be locally-Euclidean. Accordingly, 
the quantities and equations in the following sections are given locally within a small Ax or At from each point. 

The former a priori property on smoothness restricts the application of the proposed approach, while the latter a 
priori property on locality enlarges it. For example, it remains in question whether a curved surface can successfully 
model the atrium which has a substantial thickness in some parts of it, but the second a priori property supports 
the usefulness of this approach since the proposed approach can still be applied to the part of the atrium which has 
a negligible thickness or where the wave pattern is always periodic in the direction of the surface normal vector. 
Similarly, it remains in question whether we can find anisotropy satisfying the specific orientation and smoothness as 
suggested, but the second a priori property may also ensure the applicability even for the complex anisotropy of the 
atrium. 



C. Notations and definitions 

In the rest of the paper, the (double) subscript and (double) superscript index means that the corresponding 
quantity is a covariant and contravariant tensor, respectively. In principle, scalars are represented with the regular 
font without any index and vectors being represented with the bold font, but the scalar variable with the upper index, 
as the component of vector, is also used to indicate the corresponding quantity is a vector. The lists of frequently 
used notations and definitions are displayed in Table |l] and |ll| 





TABLE I: Definitions and Notations I 




ith trajectory, jth point in the ith trajectory 




jth wavefront, ith point in the jth wavefront 


X 


Afiine parameter to parameterize the trajectory 


n 


Selector parameter to parameterize the wavefront 


yi 


Trajectory with sufficiently small curvature 


y2 


Wavefront with sufficiently small curvature 




Orthogonal curved axis defined on the curved element 



The remaining sections are organised as follows: Section 2 defines the trajectory and the wavefront in cardiac 
excitation propagation. Section 3 describes the FitzHugh-Nagumo equations on curved surfaces and, in Section 4, 
using the properties of cardiac excitation propagation and a proposed supposition, we derive the Eikonal equation on 
anisotropic curved surfaces. The motivations of the relative acceleration approach are explained intuitively in Section 
5 and the hypothesis follows to yield the relative acceleration equations in Section 6. Section 7 displays some examples 
of conduction failure which are coincident with the analysis from the relative acceleration approach and the validation 
by some previous reports and computational modelling. For a practical application. Section 8 displays the analysis 
and computational modelling of cardiac excitation propagation around the pulmonary veins with various anisotropy. 
Comments and future works are provided in Section 9. 

II. TRAJECTORY AND WAVEFRONT IN CARDIAC EXCITATION 

In cardiac electrophysiology, the wavefront is well-defined as the isochrone of the depolarisation phase in cardiac 
action potential. The direction of the wavefront is consequently along the wavefront. On the other hand, the 
propagational direction is in the direction of the excitation propagation. In view of the wavefront, the propagational 
direction at a point p, which is located in the wavefront, represents the infinitesimal change of p during an infinitesimal 
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FIG. 1: The smooth curve Pi is the trajectory with an affine parameter A and the smooth curve Sj is the wavefront with the 
selector parameter n (left), yi: direction of the propagation, y2: direction of the wavefront (right). 



time interval At. Consequently, the propagational direction is orthogonal to the direction of the wavefront in isotropic 
media. For example, the velocity vector of the excitation propagation is in the propagational direction. 

Nevertheless, the concept of trajectory remains in question since the trajectory implicitly indicates that something 
is moving toward the propagational direction in the excitation propagation. Leaving what is the moving object in the 
direction of the excitation propagation for future studies, we discuss how to construct and define the trajectory in 
cardiac excitation propagation in terms of the wavefront and the propagational direction. 

To obtain the trajectory at each time and space, we first consider the wavefront Sj at time tj. Let n be the 
parameter of Sj such that Sj = Sj{n). Let's call n the selector parameter. Let = So{n) be the wavefront at time 
to. At ti = to -\- At, we will have another wavefront Si = Si{n) which is also parameterized by the same selector 
parameter n. See the left plot of Figure [ij Suppose the wavefronts 5*0 and Si are both continuous and different iable. 
For an arbitrary point So{n) on the wavefront 5*0, let the propagational velocity be v^(n). Then, the trajectory will 
be defined similarly to that in classical mechanics [3 : 

Definition: For an interval / G i?, a trajectory Po{t) is a differentiable mapping Pq : / ^ Ai^ to satisfy 



dt 



= v n = lim — —. 



For example, let Po{t) be the trajectory for the selector parameter in /c discrete time steps. For the final time T, 
let At be the time interval as At = T/k. In isotropic and homogeneous media, the propagational direction is normal 
to the wavefront, but in the presence of anisotropy, it is aligned in the direction of anisotropy. For sufficiently small 
At, the trajectory at time to + At passing So{n) is 

Po{to + At) = So{n) + v°(n)At = Si{n). 

Note that the trajectory Pn meets the wavefront Sj with the same parameter n, which is why n is called the selector 
of a trajectory. By repeating this procedure for to — At, we can construct a continuous and differentiable curve in 
a neighbourhood of S^{n) such that for a sufficiently small we have 

lim \\Pk(to + nAt) Pkito + nAt)|| < 5, for all n. 

At^O 

With the afiine parameter A = at + 6, a, b G R+, or the time that is measured by the moving body's clock, the 
trajectory is represented as P^ = P/c(A) G as well as the wavefront Sj = Sj{n) G M. What remains is to prove the 
existence and uniqueness of a trajectory for each point on the curved surface M. Consider the following proposition: 

Proposition 1: Consider a curved element that is locally-Euclidean and two-dimensional manifold. For 
any point p on the wavefront Sj that is piece-wise continuous and differentiable in M^, there exists a unique and 
piece- wisely differentiable trajectory P^ in the neighbourhood of p. 

Proof: According to the above definition of the trajectory, it is sufficient to prove that, for any point 
there exists the unique vector v that is orthogonal, in the sense of the metric gij, to the tangent vector of the 
wavefront Sj. The consideration of the general metric tensor gij is due to anisotropy. Let vsj be the tangent vector 
of the wavefront Sj. The uniqueness of vsj is provided by the assumption that the wavefront is continuous and 
differentiable. Since A4^ is a two-dimensional manifold, there is a function C which maps an Euclidean element 
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ft^ G M? into the curved element such that C : fl'^ ^ M^. Without loss of generality, let q be the point in ft such 
that C{q) = p. Let Sk be the Euclidean axis of ft^ and let d/dsk = Vgk be the tangent vector of the Euclidean axis 
for = 1, 2. Then, considering the tangent vector is mapped by C such as [45 

C (vsi) C (^^^ = ^ v/e, for each k, 

we see that each dC/dsk becomes the tangent vector at p G Then, (vi, V2) constitute a linear basis of the tangent 
plane at p where the tangent vector Msj lies, by definition. Thus, by an orthogonal procedure, for example by the 
Gram-Schmidt process in the tangent plane [43^, there exists a unique v such that v is independent of wsj and 

(v,v5j) =gij. □ 

It is clear to see that the uniqueness of the trajectory in the neighbourhood of a point p is the direct consequence of 
the uniqueness of the propagational vector and the smoothness of the wavefront. For example, for an isotropic sphere 
with a point initialization from the pole, the wavefront is aligned along the azimuthal angle and the trajectory along 
the polar angle for every point. It should be pointed out that the trajectory does not have to be orthogonal to the 
wavefront everywhere because of the feasible presence of anisotropy on the surface. 

In the following sections, we derive the relative acceleration equation in order to use the trajectory for the study of 
cardiac excitation propagation. From the Eikonal equation of the FitzHugh-Nagumo (FHN) equations, we will derive 
an equation to describe the dynamics of the normal vector (defined in Section 5), which displays the changes of the 
divergence and convergence of the trajectories, on anisotropic curved surfaces. Since the derivation is mostly based 
on the diffusion operator on curved surfaces with the reaction operator which is constant along the wavefront, the 
result of the analysis is least affected by the choice of the FHN equation for the description of the dynamics of cardiac 
excitation propagation. 



III. FITZHUGH-NAGUMO EQUATIONS ON ANISOTROPIC CURVED SURFACES 

Consider a smooth curved surface M.. Suppose that is a two-dimensional manifold such that the surface can 
be divided into a finite number of regions for a one-to-one map from a simply connected region of the Euclidean 
plane. Moreover, we suppose M. is an orthonormal surface with two orthogonal axes (xi, X2). The use of orthonormal 
surfaces has been standard in studying waves on curved surfaces for simplicity [20 [32 , but we will make the problem 
more tractable by using a piecewise orthogonal curved surface rather than a curved surface with globally continuous 
orthogonal curved axes. Consider piecewise orthogonal surfaces from the tessellation of the surface M such as, 

M = {UeM^lM'nM^ =5}}, (1) 

where z, j is the index for the element and dj is the Dirac delta function. Suppose that the curved element 
is locally Euclidean [7] such that, in a sufficiently small neighbourhood of any point p G the surface can be 
represented in a small domain of the Euclidean space with the same metric ds"^ = gijXiXj for a given metric tensor gij 
and the curved axis Xj. Note that each M^, with the orthogonal axis x^, constitutes the surface M such that M has 
two orthogonal axes everywhere, but can be discontinuous across the interfaces of the curved elements. For example, 
the above conditions are well satisfied for the surface of revolution such as cylinder, sphere, and torus. Then, cardiac 
excitation propagation on anisotropic curved surfaces can be modeled by the following Fitzhugh-Nagumo (FHN) 
equations [T4] [T8] [38]: in M\ 

§ = Gin.vj, (3) 

where g is the determinant of the metric tensor g^j. The variable u is an activator and represents the membrane 
potential, while the variable v is an inhibitor and represents the ion channel openness or the refractoriness. The 
following analysis is independent of choice of the function F{u^ v) and G{u^ v) which models the shape of the cardiac 
action potential. Nevertheless, for computational simulations, we use the following functions such as ^36j 

F{u,v) = ciu{u — a){l — u) — C2V, a, ci, C2 G M, (4) 
G{u,v) = bi{u-b2v). hi, 62 eM. (5) 
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The constants, a, 61, 62, ci, C2, determine the excitability of myocardial tissues, but in the rest of this paper, we 
only consider these variables remain constant in A4^ and consequently in A4. Thus, inhomogeneity is not considered 
in this paper. As for anisotropy, we only consider the simple anisotropy that is aligned along one of two axes so that 
the diffusivity tensor d^", written in the orthogonal curved axis Xq/ , IS expressed as 

^/3a _ ^/3a^^^ ^^^^^ ^ I 1 if a = /3 (6) 

otherwise 



Mathematically, as shown in equation ( |T2| ), the above equality can be made possible in the orthogonal curved axis 
since the diffusivity tensor d^^ is regarded as a covariant tensor of rank 2, the same as the conjugate metric tensor 
g(3a g2] [45] [43. Note that the above anisotropy can be easily constructed since each curved element is assumed 
to be locally Euclidean. Moreover, the above property is only enforced in each element A^^, thus is not a strong 
restriction with the flexible choice of curved axis that may not be continuous across the interfaces of Ai. On this type 
of anisotropic curved surface, equation ([2| boils down to 

^ ' ' ^^gd--p-)+Fiu,v), in M^. (7) 



IV. EIKONAL EQUATION ON ANISOTROPIC CURVED SURFACES 



The FHN equation derived in equation ([7|), along with equation ([3]), models the cardiac excitation propagation 
on anisotropic curved surfaces. However, the dynamics involving motion in the direction of propagation is our only 
concern. To derive the relevant equation for the dynamics in the direction of propagation, known as the Eikonal 
equation^ we use the unique property of cardiac excitation propagation as a traveling wave solution [19] [26]. For 
instance, the membrane potential u can be expressed with a function and a wave-speed c = c{yi^y2) such as 

u{yi,y2,t') = ^l){yi - c(?/i, ^2)^', ^2), (8) 

where yi is the path of the propagation and y2 is the isochrone where the variable u is constant with respect to y2 
such as du/dy2 = 0. See the right plot of Figure [l] Also, f indicates the elapsed time from the wavefront with respect 
to the axis yi and ^2, thus t' = indicates the transition layer or just the wavefront for each A. Considering the 
definition of the trajectory and the wavefront in the previous section, we see that yi is in fact the trajectory and y2 
is the wavefront. This can be confirmed by differentiating equation ([s]) with respect to y2 such as 

dyi dy2 dy2 

which shows that di\)ldy2 = at the moving frame where = 0. The next step we are going to take is to align the 
trajectory and the wavefront with two curved axes in the element M.^ . For this purpose, the following supposition is 
proposed: 



Supposition: The curvature of the trajectory and the wavefront is sufficiently small so that the axis y^ can 
be aligned with them. In mathematics, for each there exist an orthogonal curved axis x/c G M.^ such that at every 
points of A^^, for every index j and its corresponding index /c, 

< £, for a sufficiently small £ > 0. (9) 



dyi \dxkj 



In other words, the above supposition implies that the following analysis is only valid when the wavefront is slightly 
curved. But, our analyse still remains valid to indicate the initiation of a large front curvature that leads to conduction 
failure, since the critical component in the relative acceleration equation is the acceleration term as a second derivative, 
not the velocity term as a first derivative. In fact, this supposition has been frequently used in the kinematics approach 
[14j ^llj [12j to result in omitting the high order correction terms, thus the normal propagation velocity V of a curved 
wavefront is represented by a linear function V = Vq — Dk^ where D is the diffusivity constant and k is the curvature 
of the wavefront. 

In equation ([7|), the use of the traveling wave assumption ([8| yields the following equality at the transition layer 
where f = 0: 



dip 
dyi 



1 

71 



d 

dyi 



11 

dyi 



d 

dy2 



-21^ 
dyi 
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where the hat notation is used to indicate that the corresponding quantity is expressed with respect to the generahy 
non-orthogonal yj axis. For example, the diffusivity tensor df^^ is aligned along the yj axis where dP'^ is not zero 
because anisotropy can be aligned at an oblique angle with the yj axis. This is similar for the determinant of the 
metric tensor g. After sorting the terms in the above equation, we obtain the Eikonal equation of the FHN equation 
for anisotropic curved surfaces: 



d 



11 



d{yiy 



1 d^g 
9yi 



11 



dd 



11 



dyi \fg dy2 




0. 



(10) 



Since the axis yj depends on the shape of the wavefront, the axis yj can be regarded as a moving axis for t' = and 
consequently the above Eikonal equation can be regarded that it is written in a moving axis on the curved surface. 
Also, we observe that this Eikonal equation is similar to the three dimensional Eikonal curvature equation by Keener 
but there is one difference. Due to the difference between the orthogonal curved axis and the Euclidean axis, the 
component {1 / \/g){d^/g / dyj)d^^ ^ j = 1, 2 in the above equation replaces the following term in the three dimensional 
Eikonal curvature equation: 



dyi dyj \dxk 



d^\ 



J 



1, 2. 



Equation (10) is well defined in A^, but is not convenient for the further analysis because the axis yj is neither 
orthogonal or fixed. The axis yj depends on the behaviour of the propagation such as the initialization of the 
propagation and subsequently the shape of the wavefront and the time variable. Thus, it is inconvenient to express 
all the tensors with the axis yj. Instead, we express every tensors on a fixed and orthogonal curved coordinate axis 
Xa which is independent of the direction of propagation or the wavefront. 

Consider that the diffusivity tensor d^^ is expanded with respect to the curved axis Xk such as [42] [45] 



^-^ ^-^ dXk dx^ ^^r^- f)nri.'' 



/e=l j = l 



fc=l 



dxi, dxk ' 



(11) 



where the second equality is obtained from our choice of anisotropy that is only aligned along one of the curved 
axes (equation ([6|). Moreover, consider that the diffusivity coefficient d^^ can be also transformed by the Euclidean 
axis x"^ . Without a loss of generality, we can let the Euclidean axis x^ be constructed such that the anisotropy 
in the Euclidean space with the axis is aligned along one axis only. Note that this is possible because of the 
aforementioned properties of the surface that is locally Euclidean. Then, we obtain 



E 

k=\ 



^l^jjdxkdxk 
^ V dx^ dx^ 



dyadyi ^ y^^kgkk dyg dy(3 
dxk dxk f-f dxk dxk ' 

k—l 



(12) 



where we used 



d^^ = d^ det 



dx^ 
dxk 



and g^^ = Yl 



dxk dxk 
dx^ dx^ ' 



For example, let's introduce a new variable Ak and a new tensor defined such as for /c = 1, 2, 



dyi 

dxh 



and l]^=d^g^^Ak. 



Then, we can express d^^ and d'^^ in equation (10) as follows: 

2 



d^^ = ^ A^ and d^^ = ^ 



k ( dy2 

U 

k=l /c=l 

As for ^/g^ the transformation rule by multiplying the determinant of the Jacobian applies as [28] 

dxn 



= where J = det 



(13) 



(14) 



(15) 
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As a result, by substituting the equations (14) and ( |T5| into equation ( [To) ), we obtain 

d{yir [Vy dxi ^ dx2 } dyi 



(16) 



where the speed variable c is newly defined as 



k=l 



J Sx/e dyi\dxk) dy2\dxk) 



}■ 



Moreover, referring again the supposition (|9| to imply that, for the moving axis y^^ we can find a fixed curved axis 
Xk which is not a large deviation from the orthogonalized axis yk^ we observe that the following components are 
sufficiently small such as, for a sufficiently small (5 > 0, 



1 dJ 

J dxk 



< S, and 



d 






\OXkJ 



< (5, for k, j = 1, 2. 



Consequently, the speed function c can be approximated as c, independent of the geometrical factors such as d^, A^, 
and i.e., c ~ c. 



V. HYPOTHESIS OF THE RELATIVE ACCELERATION FOR CONDUCTION FAILURE 

Before proceeding further to the relative acceleration equation, in this section, we will briefiy describe the meaning 
of the relative acceleration in terms of the trajectory and explain how the relative acceleration can be translated 
into the conduction failure of cardiac excitation propagation. The goal of this section provides the backgrounds and 
justifications of the proposed hypothesis, which plays a critical role for the rest of the analysis, from discrete models 
similar to cellular automata [9 . Based on the crucial characteristics of the excitation propagation, the validity of the 
hypothesis is obvious independent of the scale of the system, but the proof of it remains challenging and beyond the 
scope of this paper. Thus, we just leave it as a hypothesis for later rigorous validation. 



A. Meaning of the relative acceleration 

Let Po be a trajectory as previously defined. In the neighbourhood of Pq, consider the other trajectories Pk for 
k = ±1, 2, 3... (Figure[2|. All the trajectories are given the same affine parameter A so that Po(A) and Pfe(A) lie on the 
same wavefront Sx. In other words, the relative acceleration of Pq measures how Po(A) advances with respect to other 
P/c(A). Two kinds of relative accelerations can be considered: the first relative acceleration is in the propagational 
direction and the second relative acceleration is in the direction of wavefront, or in the normal direction. Relative 
acceleration in the propagational direction is obvious. Suppose that the trajectories are all parallel. If the rate of 
changes of ||Po(A) — Po(A — 1)|| is more than the first order compared to the other trajectories, then Pq is said to be 
relatively accelerated in the propagational direction. 

On the other hand, the relative acceleration in the normal direction should be considered for three different cases: 
suppose that the rate of change ||P/c(A) — P/c(A — 1)|| is the same for all the trajectories. First, when all trajectories are 
parallel as shown in Figure [2]A, the change of the separation vectors (as defined below) along each trajectory is zero, 
thus there is no relative acceleration in the normal direction. Secondly, if the trajectory Pq linearly diverges from 
the other trajectories as shown in Figure [2^, then the first derivative with respect to A is nonzero, but the second 
derivative is still zero. Consequently, there is no relative acceleration in the normal direction, either. However, if the 
trajectory Pq diverges quadratically from the other trajectories as shown in Figure |2p, none of the first or the second 
derivatives is nontrivial, thus the trajectory Pq is said to be relatively accelerated in the normal direction. 

More rigorously in mathematics, we use the following definition of the relative acceleration which is equivalent to 
what has been explained: 

Definition: Let n^ G be the separation vector to indicate how two trajectories are separated for the 
same A of the trajectories (excuse the abuse of the letter n (selector parameter) again). If the separation vector n is 
defined such as 



n{i,k)= lim Sk{i ^ An) - Sk{i) ^ 

An^O 



(17) 
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then the relative acceleration [30^ can be expressed as, with n = (n-^,n^,n^), 




FIG. 2: Three different cases of the relative acceleration in the normal direction, is the separation vector connecting the 
same A for the trajectories. 



B. Insights from discrete models 

To qualitatively explain how a large relative acceleration can be interpreted as conduction failure for the justifi- 
cation of the following hypothesis, we first consider stopping conditions in a discrete model consisting of individual 
myocardial units that represent cells or tissues. Motivated by the widely known facts on myocardial cells or tissues, 
we give the following properties to myocardial units: 

(i) One unit faces multiple units such as hexagonal packing. 

(ii) One unit can excite a limited number (A/'^ax) of neighbouring units and when they try to excite more than A^max, 
none of the neighbouring units reach their threshold potential (TP) for excitation. 

(iii) After a unit is excited, it undergoes the recovery process, called the refractory period (RP), and is not excited 
during RP. 

For the structural property (i) of myocardial tissues as reported in ref. [24] [35], the property (ii) provides 
the unique characteristics of excitation propagation to restrict the maximum number of excited cells that one cell can 
excite, known as impedance mismatch or sink-source mismatch [37] [46] . Property (iii) is also well known. Depending 
on ionic currents of (7a^+ and the minimum time for the recovery of excitability is expressed as a function of 
RP that is also correlated with action potential duration (APD) [33 . 

Figure shows myocardial tissues modeled as hexagonal packing. The geometric-free discrete representation is 
shown in Figure [3^, where a dark colour indicates the cell is excited and white colour indicates the cell is in excitable 
state. Figure |3p shows the wavefront to represent the line of excited tissues at each time. 

For the hexagonal packing of myocardial units as shown in Figure [3]A., we can also represent it with a geometric- free 
discrete model as shown in Figure [3^ in which the number of units in each layer reflects the shape of the hexagonal 
packing. Starting from (A1,A2) units, the excitation propagates to (B1,B2,B3) units, (C1,C2,C3,C4) units and so 
forth. Note that the first letter of the unit indicates that the same layer must be excited in all of the units and the 
excitation occurs only in the alphabetical order, i.e., a C-unit cell cannot excite a B unit or a C unit. The goal of this 
discrete model is to expand the concept of the unit from cell to tissue, because when a cell satisfies the aforementioned 
properties, a tissue being a group of cells also satisfies the same properties if we neglect the time difference of excitation 
within a tissue. 

No matter what this one unit in the discrete model represents, if we connect excited units at every time interval, 
we obtain the wavefront of cardiac excitation propagation. This wavefront represents a series of myocardial units 
that are in the process of depolarisation. Consequently, the trajectory can be defined as if the cardiac wave is a 
physical wave, that is, a collection of moving particles. Upon these representations and notions, we describe the 
stopping conditions of the excitation propagation as the following: Without a loss of generality. 

Let three be the maximum number of units one can excite in the model of hexagonal packing^ i.e., A^max = 3. 
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The first stopping condition is obtained directly from property (ii). In Figure [Sp, if the Al and A2 units at- 
tempt to excite the B-units, then each unit needs to excite an average of 3.5 units. Consequently, the electric 
potential of all B-units reaches below the TP and, consequently, fails to be excited. In brief, the excitation 
propagation stops or conduction failure occurs if: 

The average number of units to be excited exceeds the maximum number of excitable units. (SC-i) 

This stopping condition has often been cited when excitation propagation faces abrupt tissue expansion, such 
as a narrow cell strand connected to a large rectangular cell [37]. In the same context, we can explain why a small 
gap between ablation lesions, burned and dead myocardial tissues by catheter, does not allow excitation leakage. It 
is important to note that the excessive number of the B-units represents the geometric distribution of myocardial 
units, which can be interpreted as geometry in physics. 

On the contrary, less attention has been given to the second condition. In Figure [3^, the geometry of hexagonal 
packing is the same as that of normal propagation in Figure [3^, although the Al unit propagates rapidly only along 
the Bl unit, and finally to the CI unit. This type of propagation is common in myocardial fibers or myocardial sheets, 
which we generally call anisotropy [47 . Because the CI unit is the only excited cell in C-units, the CI unit attempts 
to excite all the five D-units. Let Tl be the time when the CI unit uses all of the electric potential for the D-units. At 
the same time, the A2 unit proceeds to excite the (B2, B3) units and subsequently the (C2, C3, C4) units. Let T2 be 
the time when the (C2, C3, C4) units are all excited. There is a critical difference concerning the times Tl and T2: if 
T2 is earlier than Tl, that is, if the (C2, C3, C4) units also can contribute to exciting D-units along with the CI unit, 
then all the D-units can be excited and the propagation continues. However, if Tl is earlier than T2, that is, if when 
the (C2, C3, C4) units are excited, D-units are already in RP so that the (C2, C3, C4) units have no units to excite 
and propagation stops at the D-units. In other words, excitation propagation also stops, or conduction failure occurs if 

The excitation time of other C-units exceeds the propagation time from the CI unit to the D-units consisting 
of more than the maximum number of units one unit can excite. (SC-ii) 

To simplify the above condition, we introduce another property of myocardial units as: 

(iv) The excitation time by a myocardial unit is roughly proportional to the number of cells that one cell has to excite. 

With the property (iv), the time Tl only depends on the propagational time from the Al unit to the CI 
unit, because the time required to propagate from the CI unit to the D-units is fixed. Moreover, the ratio of velocity 
is at maximum four times higher on fibers than the normal excitable media Thus, it is reasonable to say that 
the time Tl is relatively fixed. On the other hand, the excitation time of other C-units depends on the geometry 
of the B-units and the C-units other than the Bl unit and the CI unit. Thus, the time T2 varies widely and is 
dependent upon the geometry of the C-units and the D-units. 

These discrete models can be naturally interpreted in continuous wave forms. This procedure can be performed 
by drawing a line for excited units at each time which represents the front of the excitation propagation, called a 
wavefront in light of the waves. For example. Figure [sjC displays the wavefront (solid line) and trajectories (dotted 
line) of the excitation propagation for the discrete model of Figure [3^. In waveforms, the shape of the wavefront and 
trajectories now refiect the geometry of myocardial tissues. 

Similarly, the wavefront and trajectory of the first stopping condition is displayed in Figure [3^. A high ratio 
between excited cells and excitable cells leads to abrupt expansion of the distances between the trajectories. In 
classical mechanics, this phenomenon can be described as a large relative acceleration in the normal direction, that 
is, a non-uniform rate of change of speed along the wavefront. In homogeneous media, this means that the distance 
between particles continues to diverge along the wavefront as the wave propagates. On the other hand. Figure [3p 
represents the wavefronts and trajectories of the second stopping condition, but instead of the spreading-out of the 
trajectories, we observe a relative acceleration in the propagational direction. 

In other words, without regarding the directions of relative acceleration, conduction failure (stopping condition) 
for cardiac excitation propagation can be expressed as the following hypothesis: 

Hypothesis: A sufficiently large relative acceleration of cardiac excitation propagation causes conduction 
failure in the homogeneous and anisotropic myocardial tissues. 

Note that this hypothesis may not hold in inhomogeneous media such as damaged or malfunctioning myocar- 
dial tissue, but has no restriction on the dimension of space or the changes of the conducting velocity. 
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FIG. 3: Representation of stopping conditions for conduction failure in individual myocardial cells and in continuous waves. 



TABLE II: Definitions and Notations II 



J 



Jacobian from the ^/j-axis to the x^-axis 

the component of the propagational vector in Xk axis, or dyi/dxk 
tensor that is equivalent to d^g^^Ak 

ELi(yV9)d{V9^')/dxk 



VI. RELATIVE ACCELERATION EQUATION 



With the proposed hypothesis, what remains is to obtain the critical factors for the relative acceleration in the 
excitation propagation. We begin with regarding the ?/i-axis as the trajectory in the ambient space of or R^. 
This procedure can be considered as the projection of yi into the ambient space of a higher dimension [8] [30] . In the 



Eikonal equation (10) which is defined in A^^, we extend the curved surface M.^ in the neighbourhood of each point 
in the following way: consider a curved surface as a two dimensional submanifold embedded in three dimensional 
Euclidean space x\ The function ^^(^1,^2) can be extended to a function in a tubular neighbourhood 

of the surface where the directional derivative of along the surface normal is zero. Then, is defined on an open 
subset around A4^. More rigorously, let JV be the parameter in the direction of the surface normal that is also unique 
due to the assumptions on A^^. Then, we can expand the function tj; outside the surface such that dt/j/dJV = in the 
neighbourhood of p G A4^. Let Yi : ^ be the coordinate map for the ?/i-axis in the ambient space R^. With 
the affine parameter A and the selector parameter n, we introduce the following notation as, 

P„(A) = Yi(A,n) = (7^, 7^, 7^)eK3, 

where 7^ is the component of Yi for the Euclidean axis x*. Moreover, since '0(xi,X2) was extended to x^, x^), 

we can express the differentiation with respect to yi by the chain rule such as 
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Or, using that dx^ /dyi is the ith component of the tangential vector of the path Yi that is the same as 97^/9 A, we 
obtain 



dyi 



E 



i=l 



dx^ dX ' 



(19) 



Moreover, differentiating again with respect to yi yields 

^2 / 3 



dyidy 



= E 



dx^ 



E 



dx^dx^ dX dX 



(20) 



By substituting the equalities ( |19[ ) and ( |2Q[ ) into equation ( |16| ), we obtain 

1 di^V^) , 



KM 



H 
dX 



V J 



dx'^dx^ dX dX 



0. 



where the summation notation is used for easier reading where the index i and j is summed up to three and the 
index k is summed up to two. To obtain the relative acceleration equation, we differentiate the above equation with 
respect to the selector parameter n along the wavefront. Considering that the reaction function F is constant along 
the wavefront, leading to dF/dn = 0, we derive the following equation with summation notation: 



dx^ 



d^n' 5(AfeU^) 



dX^ 



dn 



'dX 



dx^dx^ 



d_ 
dn 



1 dj^V^ 
dxk 



}- 



(21) 



where we used the variable n* for the separation vector as defined in equation (17) and v'^ for the tangent vector of 
the trajectory Yi, which are defined as follows: for 1 < i < 3, 



dn 



and 



dX 



Each upper index indicates the quantity is the component for the Cartesian coordinate x*, for example, n = 
(n^, n^, n^). In addition, we used the interchangeability of the differentiation [30] in the direction of n and A 
for the above equation such that 



d_ 
dn 



d_ f d_ 
dn \ dX 



d_ 
dX 



In 



d_ 
dX 



d_ 
dX 



dn 



d^n' 
~d>?' 



In equation (21), let us consider the case only when the quantity inside bracket is zero, independent of dijj/dx^. 
Without a loss of generality, we obtain for each index i. 



~dX^ 



1 



dn 



, dv'^ 
dn 



dE fdv' d^ip 

— \ — _ ■ _ ■ ^ 

dn V dX dx^dx^ 



d{cv'^) 
dn 



(22) 



where we introduced the new parameter E and G that are defined such as 

2 2 



k=l 



k=l 



dxk ^ ^_^dxk^^ 



Note that each component in the right hand side is only dependent on the coordinate axis, not the initiation of 
the propagation or the the time variable. For the above components for the relative acceleration equation (22), 
we notice that the last two components are relatively small compared to other components with the following 
reasons: first consider the following term (d'^tjj/dx'^dx^). Since tjj is the membrane potential, this component, 
differentiated twice with respect to Xj, represents the static electric force produced by the shape of cardiac action 
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potential. Considering the smooth and uniform shape of cardiac action potential and the relatively weak amount 
of the electric voltage, the static force is relatively small compared to other dynamic forces for the excitation 
propagation. Moreover, considering the supposition and the fact that cardiac excitation propagation is slowly 



varying as stated in ref. [25 , we see the fourth component of equation (22), d{cv'^)/dn, is independent of E and G 



and also relatively small. As a consequence, we propose the following proposition with the aforementioned hypothesis. 
Proposition 2: Suppose that the curved surface M is homogeneous and anisotropic. If the magnitude of 



equation (22) is sufficiently large at a point p G Al, then the excitation propagation stops at p leading to conduction 



failure in the neighbourhood of the point p. 



Proof: With the hypothesis in Section 5, the relative acceleration equation (22) is obtained, as shown in 



Section 6, from the Eikonal equation (16) of the FHN equation on curved surfaces. □ 



A large magnitude of equation (22) may also mean a large relative deceleration, but the use of the anisotropy that 
is larger than 1.0, which means the faster conducting velocity than the normal tissue, rules out a large deceleration. 
Note that the relative acceleration does not have to occur in a large area simultaneously for conduction failure, but 
it suffices to occur only at a point since the propagational failure can diffuse subsequently from one point to other 
points along the wavefront. This can be easily pictured by the game as explained in Figure |3] 

VII. EXAMPLES 

In this section, in order to demonstrate the validity of the proposed relative acceleration approach, some conditions 



for conduction failure from the relative acceleration equation ( 22 ) are validated by simple computational simulations 
on curved surfaces and/or by the previously known results from the kinematics approach. 

For the computational simulations, Nektar++ |40], a C++-object-oriented partial differential equation solver, is 
used with partial modifications according to the validated scheme, known as the method of moving frames [8 , to 
solve partial differential equations on anisotropic curved surfaces [10]. Consider the mixed formulation of the FHN 
equations such as 

Ou 

^ = V • q + F(ii,v), q = dVii qeM, 

^ = G{u,v), deM, 

where F{u^v) = ciu{u — a)(l — u) — C2V, G{u^v) = bi{u — 62'^), a = 0.13, bi = 0.013, 62 = 1-0, ci = 0.26, C2 = 0.1. 
For the above FHN equations, 1 time unit is equivalent to 0.63 ms and 1 space unit to 0.99 mm [36 . Note that both 
of q and the diffusivity tensor d lie in A4. Integrating the above equation with respect to test functions, we obtain 
a weak form of the above equations which is expressed with the orthogonal moving frames [7] of the unit length 
such as 

/ qmie"^ ■ V(p)dx - / qmi^"^ ■ n)(pds + / F{u, v)dx, 
I qm{dr)~^^dx = - [ (V • {(pe'^))udx + / (e^ • n)ipuds, 





du ^ 


1 




/ 










dv , 


f 


-V>dx 



I 



G{u, v)dx, 



where d^ = d - and the tilde sign represents the approximated value of the corresponding quantity, or known as 
flux in the context of the discontinuous Galerkin methods. 

In the following examples on various curved surfaces, we will display how mathematical analysis from the relative 
acceleration equation coincides with the computational results and how the curvature of the geometry changes the 
impact of anisotropy for conduction failure. For easier reading, each example is organized such that Problem means 
the setting of cardiac excitation propagation on a curved surface, RA-analysis means the analysis and predictions 



from the relative acceleration equation (22), Computational modelling and validation means the validation by 



computation simulations in the above computational scheme and by some references, respectively. 
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A. Anisotropic plane 

Problem 1: Consider a block of anisotropy in the middle of the plane. For the plane with —100 < x < 100 
and —100 < y < 100, let the block be located in —50 < x < 50 and —50 < ^ < 50. Let the anisotropy be aligned 
with the direction of the x-axis, denoted as x -anisotropy which expresses the diffusivity tensor d as d = (i^x 
in the block of anisotropy and = 1.0 everywhere. See Figure [4]A.. Also, let the excitation propagate in form 
of a plane wave in the —x direction front the right wall and, consequently, the wavefront is in the direction of the ?/-axis. 

RA-analysis: Since these conditions mean that 

^ — 

ox 



since g"^"^ = = I and A^^ = 1, Ay = 0, 



the relative acceleration equation (22) boils down to 



1 



dy dX dxdy 



d{logd^) dv^ 



1 a^d^ 



dy dX dxdy 



(23) 



where we used n = y and dv'^/dn = 0. From the above equation, it is obvious to see that a large relative acceleration 
can be achieved only by increasing the magnitude of anisotropy in the plane, especially at the interfaces of anisotropy 
where d'^d^ /{dxdy) is large. This result is compatible with the break-up conditions drawn from the kinematics 
approach as shown by Morozov et. al. [31j . 

Computational modelling: In the computational simulation, we similarly observe that a large relative ac- 
celeration can only occur at the interfaces of the anisotropy block, i.e., in the line of y = ±50 for the first component 
and two points (50, 50) and (50, —50) for the second component. See Figure [sjA. for the distribution of d^ and 
d{\ogd^)/dy along the ?/-axis. With d^ = 4.0, conduction failure does not happen as shown in Figure [4^ and[4]C. 
However, increasing the magnitude of anisotropy up to d^ = 10, though this large magnitude seems to be unrealistic 
in the biological system, leads to conduction failure in the block of anisotropy because the relative acceleration has 
increased significantly at the interface. Note that if the anisotropy changes slowly such that d'^d^ /{dxdy) remains 
bounded, then no conduction failure would be induced by the anisotropy. 



Problem 2, RA-analysis and validation: Conduction failure can also occur even inside an anisotropic 
block. If the direction of anisotropy is oblique to the propagational direction, then conduction failure can happen 
inside anisotropy. Let d be the magnitude of anisotropy and be the angle with respect to the direction of the 
excitation propagation. For the same planar propagation as before, the relative acceleration equation is expressed as 



d{\og d) 
dy 



dO] dv' 
dtanO—- — - 
oyj oX 



1 d^d 
_d dxdy 



tan^ 



dlogdde dlogddO 
dx dy dx dx 



dO dO 
dx dy 



tan^ 



d^O 
dxdy 



which implies that the excitation propagation can also be blocked by significant variations of the direction of anisotropy. 
For example, let the magnitude of anisotropy {d) be constant in the domain. Then, the above equality reduces to the 
following equation: 



-dtdbiiO 



dO dv' 
dy dX 



dO dO 
dx dy 



d^e 

dxdy 



(24) 



The above equation confirms that the relative acceleration can also be significantly increased by the variable angle 
of anisotropy as well as the magnitude of anisotropy. This result coincides with Davydov et. al. [17 which displayed 
that the chiral anisotropy, which is aligned along the circumferential direction of circles, can break up the excitation 
propagation. But, the above equality suggests a host of possibilities for conduction failure by variations of the direction 
of anisotropy. 



B. Anisotropic sphere 



Problem 3: The relative acceleration equation (22) can be similarly applied to the propagation on anisotropic 
sphere. But, there is a difference. Contrary to the propagation in the plane, the planar propagation is impossible. 
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FIG. 4: For — 4.0. Anisotropy block is located in the middle (A). After initiating from the right wall, the membrane 
potential {u) at T = 300.0 (B), T = 500.0 (C) from the right waU. 




FIG. 5: For (T = 10.0. The distribution of (F (solid line) and d{\og(r)/dy (dashed hne) along the ^/-axis (A). After initiating 
from the right wah, the membrane potential {u) at T = 200.0 (B), T = 300.0 (C). 



Instead, the propagation should be point-wisely initiated, for example initiation from the north pole. Geometry tells 
us that the point-initialized excitation follows the line of the polar angle {0) and the wavefront is in line with the 
azimuthal angle (^). For sphere of radius r = 50.0, consider a block of anisotropy located in — 7r/4 < ^ < 7r/4 and 
—0.6 < < 0.6. See Figure [6|^. Let the direction of anisotropy be both along the 6- and (j) axis in general, but 
because of the propagational direction, only 6>-anisotropy will contribute to the changes of the propagation. 



RA-analysis: Geometrically, these conditions imply that 



E ■ 



^— - + ^cot6', smce g 



1 



and A^ = 1, A^^ = 0, 



thus, the relative acceleration equation \22\ is expressed as 

d^n' _ d{\ogd^) dv' ( 1 d^d^ d{\ogd^) 



dX^ d(t) dX \d^ 

where we used, as proved in Appendix A, for any /c, 

1 



cot e 



(25) 



The relative acceleration equation (25) is similar to the equation (23) for the plane, but the following additional 
component appears: 



d{\og(f 



cot 0. 



(26) 



The above component represents that the relative acceleration depends on the location of anisotropy with respect 
to the point of initialization. For example, for the point of initialization at the pole, let be the azimuthal angle 
of the sphere from the pole. Then, the magnitude and sign of the above component changes as the location of 
anisotropy changes. Figure [tJ^ displays that the above component changes signs at = it/ 2 and the magnitude 
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of it increases as gets closer to or tt. Consequently, since the change of signs of the above component can 
increase or decrease the total magnitude of the relative acceleration, equation (25) implies that conduction failure 
also depends on the location of anisotropy with respect to the location of point initialization. In the perspective 
of the trajectory, this analysis is no surprise. With point initialization at the north pole, the trajectories diverge 
in the northern hemisphere up to the equator and from the equator the trajectories converge in the southern 
hemisphere. Thus, the use of anisotropy to increase the divergence of trajectories will add its relative acceleration 
only for the trajectories in the northern hemisphere and will decrease it for the trajectories in the southern hemisphere. 



Computational modelling: To confirm the above analysis, place a block of anisotropy such that the first 
interface that the propagation meets is located at = 7r/4 distance from the point of initialization, i.e., in the 
northern hemisphere. As shown in Figure [6^ and [6p, the block of anisotropy with a sufficiently large magnitude 
causes conduction failure. On the other hand, move the point of initialization away from the block of anisotropy such 
that the first interface that the propagation meets is located at ^ = 37r/4 distance from the point of initialization, 
i.e., in the southern hemisphere. Then, the block of the same anisotropy, with the same magnitude of anisotropy as 
before, does not cause conduction failure as shown in Figure [7^ and[7p. This phenomenon seems to be well explained 
by the change of signs of the component (26) representing the convergence and divergence of the trajectories. 




FIG. 6: For = 10.0. The block of anisotropy is located in the middle 7r/4 <0 < 37r/4 on ; 
from the north pole, the membrane potential {u) at T = 300.0 (B), T = 400.0 (C). 



sphere (A). After point-initialized 




FIG. 7: For d = 10.0. Magnitude of the relative acceleration component as shown in equation (26) with respect to (A). The 
block of anisotropy is located in the bottom Stv/A < 6 < tv. After point-initialized from the north pole, the membrane potential 
{u) at T = 600.0 (B), T = 700.0 (C). 



C. Anisotropic torus 

Problem 4: The excitation propagation in torus is strikingly different from the propagation in the plane or sphere 
in the sense that we can not find one axis to represent the direction of propagation. For example, reminding that 
Ak = dyi/dxk, Aa^ = 1 and A^ = for the planar propagation in the plane and A^* = 1 and A^ = for sphere with 
point-initialization. However, A^* and A^ for torus is not always constant. With radius of the meridian r and radius 
of the great circle let be the axis around the meridian and (/) be the axis around the great circle. Then, a torus 
can be parameterized such as 

X = {{R-\- r cos 0) sin (j)^{R-\- r cos 0) cos (/>, r cos ^), 0<^<7r, — tt < (j) < tt. 
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FIG. 8: Components in equation (27) 
(C) for r = 0.4 (solid line), r = 0.2 (dashed 



vs. 0. T] = 7r/4 and R = 1.0. dE/dn versus 6 (A), G versus 6 (B), and 
line), r = 0.1 (dotted line). 



Then, no matter where or how the excitation starts, and A^i is always neither constant or zero. This fact makes 
the analysis complicated, but to make the problem simpler, we first consider there is no anisotropy. 



RA-analysis and validation: From the relative acceleration equation, with 



E = 



A? 



A 



smce 



(i? + r cos i9)2' 



9 



and G = 



A6I sin 



OAs 



1.0, we obtain 
1 



dAe 1 ^ 

dO r{R -\- r cos 0) d(j) {R-\-rcost 



and 



(i? + r cos (9)2' 



Accordingly, dE/dn and dG/dn in the relative acceleration equation (22) are expressed as 

2r sin 



dE 


_ 1 d{Aj) 


dn 




dG 


_ 1 d^Ae 


dn 


dOdn 



•A 



dO 

{R + rcosO)^ dn^ 
sin dAe Aq R cos 9 



r{R -\- r cos 0) dn 
d^A^ 1 ^ dA^ 

dnd(j) {R -\- r cos 0)'^ d(j) 



r {R -\- r cos 0)'^ 
2r sin 



{R- 



r cos t 



To make the problem even simpler, we suppose that the wavefront is at the angle of 77 with respect to such that 
Aq = COST] and A^ = sin 7^. Let Aq and A^ be constant along the wavefront. Then, each relativity component boils 
down to 

_ sinO dE .9 2rcos0 



= — cos a 



dG 
dn 



r{R + r cos^) 
cos a R cos 6 -\-r 



7—= cos asm a—— 
dn [R 



r cosOy 



{R- 



r cos 



\2' 



(27) 



With T] = 7r/4 and R = 1.0, Figure [s] displays the magnitude of each component versus 0. Note that each component 
displays strong dependency on the meridian r for the fixed R. In other words, the above relative acceleration equation 
for isotropic torus implies that the relative acceleration can be significantly changed by resizing the meridian r. 
This result is coincident with the kinematic analysis by Davydov et. al. [14 which demonstrated that the critical 
curvature for breaking- up the wave can be generated by adjusting the ratio R/r of torus. 

Problem 5: Consider anisotropy on torus. But, for the sake of simplicity, we suppose the propagation fol- 
lows locally the (j) axis only in a small region of torus. Let R = 100.0 and r = 50.0. Let the center of torus is located 
at (0, 0, 0). The block of anisotropy is placed in the area of —50.0 < x < 50.0 and 0.0 < y < 150. If the propagation is 
point-initialized at ( — 150.0,0,0), then the propagation approximately follows the (/)-axis in the area of the anisotropy 
block. 



RA-analysis: With the ^-axis and the 0-axis as defined previously, for Aq = and A^ = 1, we obtain 

1 dd^ 



E = 



(i? + r cos i9)2' 



and G = 



(i? + r cos l9)2 
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Accordingly, the following relative acceleration equation is obtained: 



d^n' _d{\ogd^)dv' d^d^ ^ 2rsin(9 
dX'^ dn dX dndcj) R-\-r cos 6 



dO 
dn 



(28) 



where we used d(j)/dn = according to the assumption that the propagation follows the 0-axis in the area, thus n = 0. 

The first two components are the same as the relative acceleration generated by anisotropy in the plane, but the 
last component is an additional component for anisotropic torus. Because of this last component, the behaviour of 
cardiac excitation propagation displays the following unique phenomena: the first phenomenon is that the relative 
acceleration depends on the angle. Considering 2r sin 0/{R + rcos^) increases as approaches 7r/2, it can be 
predicted that there is a larger relative acceleration in the area where is closer to 7r/2. The second phenomenon is 
as follows; because the last term adds additional relative acceleration, the critical magnitude of for conduction 
failure is slightly less than that for anisotropy in the plane. This is possible because all the components in equation 
(28) are in the same sign when d^ is larger than 1.0 which is the choice of our anisotropy. This result seems to 



support the conjecture that the curvature of geometry can increase the effects of anisotropy on curved surfaces. 



Computational modelling: The predictions of these unique phenomena in torus can be confirmed in com- 
putational simulations. Consider a torus of R = 100.0 and r = 50.0 centred at (0, 0, 0). Let the block of anisotropy be 
located in —50.0 < x < 50.0 and 0.0 <y < 150.0 as shown in Figure [9|^. In Figure [9^ and|9p, when the propagation 
is approximately in the direction of the (j) axis, the anisotropy d^ = 8.0, which is less than d^ = 10.0 for the previous 
cases, generates the largest acceleration in the area where = 7r/2. But, the block of anisotropy fails to generate a 
sufficiently large relative acceleration for the area where is small. This agrees with the first predicted phenomenon 
from equation (28). However, the d'^ anisotropy cannot block the propagation in result. This could be just a matter 
of direction of the propagation. Figure 10 B and 10 3 displays that additional anisotropy in another direction d^ can 



add additional relative acceleration in the area where is small and consequently can stop the propagation. Note 
that in the plane, whatever the direction of the propagation is, anisotropy with the magnitude of 8.0 cannot block 
the propagation for conduction failures. This confirms the second predicted phenomenon by the relative acceleration 
equation (28). 
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FIG. 9: Location of the anisotropy block (darkened) with d'^ = 8.0 (A). After point initialized at the rightmost point of the 
above torus, the membrane potential (u) at T = 500.0 (B) and at T = 1000.0 (C). 




FIG. 10: Anisotropy block with — 8.0 and d^ — 8.0. After point initialized at the rightmost point of the above torus, the 
membrane potential {u) at T = 800.0 (A), T = 900.0 (B) and T = 1000.0 (C). 
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VIII. APPLICATION TO THE PULMONARY VEIN WITH ANISOTROPY 




FIG. 11: Modeling of the PV by a surface of revolution. The (/)-axis is in the circumferential direction of the column and the 
^-axis is its orthogonal to it. 

In this section, we apply the relative acceleration approach to an anisotropic surface which models a part of the 
atrium with a pulmonary vein (PV) covered with anisotropy. The uniform smoothness of the structure, the negligible 
thickness of the atrium and the continuous alignment of anisotropy may be excessively simplified from the perspective 
of realistic cardiac anatomy, but we try to demonstrate the use of the proposed approach to provide analytical 
explanations on some reported phenomena even with this simplified model for curved surfaces. 

A. Model of the atrium with the PV 

The left atrium, of which thickness ranges from 0.5 to 3.5 mm, is a thin- wall structure of average thickness of 
1.89 =b 0.48 mm [5]. The left atrium is a lot thinner than the left ventricles which is averaged as 12.0-15.0 mm for an 
adult heart [21], but, due to the variable and non- negligible thickness in some part of the atrium, it remains a question 
whether the left atrium can be well modelled as a curved surface. Even our approach may be directly applied only 
to some area of the atrium with sufficiently thin wall unless other assumptions and simplifications are made. On the 
other hand, in the pulmonary vein (PV) which is the pathway of the oxygenated blood from the lung and is covered 
with the cardiac tissues up to a few inches from the atrium, the thickness of the cardiac tissue is thinner and less 
variable. According to Hocini et. al. [23 , the cardiac tissues of the PV are thickest close to the atrium and thinner 
toward the lung, ranging from 0.5 to 0.8 mm that is relatively constant and half times less than the average thickness 
of the left atrium. Thus, we may assume that the PV can be well modelled as a surface. 

In view of geometry, the three-dimensional shape of the PV is peculiar. Even when the shape of the atrium is 
approximated as a simple smooth sphere, the PV should be constructed as an individual surface and be rooted at the 
atrium such as a smooth column on a sphere without considering stretching it. This does not mean that the intrinsic 
curvature [27 of the PV is strikingly different from that of the atrium. In other words, the shape of the PV can be 
constructed by cutting and pasting the part of the atrium or by putting it inside out, not by stretching it. Instead, 
the PV is given with the different principal curvatures such as the Gaussian curvature. Mathematically speaking, we 
say that the PV has different exterior curvature [27] from that of the atrium. The uniqueness of the PV begins with 
the different exterior curvature, but there is one more crucial factor to consider; anisotropy. 

The real atrium has multiple pulmonary veins (4 for the human atrium and 5-6 for canine atrium), not even 
mentioning the PV-like geometries such as the inferior vena cava, the superior vena cava, the coronary sinus and 
possibly the valve annuli pL^. Consequently, it is geometrically impossible to construct a continuous distribution of 
anisotropy at all interfaces between the PVs and the atrium. This problem is similar to the mathematical construction 
of a continuous and differentiable vector field on a connected curved surface with n holes, for example, 8 holes for the 
human atrium. Thus, the discontinuities of anisotropy are highly likely to appear close to the root of each PV. 

This fact has been also observed anatomically. Hocini et. al. [23] reported that the anisotropy of the canine PVs 
abruptly changes its direction or in the mixed orientations at the root of the PVs and is aligned along the longitudinal 
direction as the vein is closer to the lung. Also, Verheule et. al. [44 observed the circumferential orientations of 
anisotropy in a vein of the canine atrium that is sectioned transversely at the various tip of the PV sleeves. Note 
that the circumferential orientation of anisotropy in the PV also presupposes the discontinuity of anisotropy at some 
roots of the multiple PVs. It still remains a question on how to construct or approximate the mixed orientation of 
anisotropy, but for simplicity we only consider the longitudinal anisotropy or/and the circumferential anisotropy. 
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In this section, we use again the relative acceleration approach to reveal the role of the peculiar geometry of the PV 
with a specifically oriented anisotropy. For this purpose, a simple surface of revolution will be used as the geometric 



model of the PV as shown in Figure 11 Consider a circular arc with radius R which is in the distance of r from 
the center of the column. The ^-axis is aligned along the circular arc and the (/)-axis is aligned along the revolving 
direction. Note that ^-axis is in the longitudinal direction and the ^-axis is in the circumferential direction. Revolve 
the circular arc around the x-axis, i.e. perpendicular to the plane, generates a smooth column on the plane which has 
the following parameterization of the surface x: For a i?, r G and < < — 7r<^<7r, 

x= (i?(l - cos l9),(i?(l - sin l9) +r) sin ^,(i?(l - sin l9) +r) cos ^), (29) 

For the metric tensors and the Christoffel symbols of this surface, refer to Appendix B. As for anisotropy, we generally 
consider the anisotropy aligned in both directions such that the diffusivity tensor d is expressed as 

d = + where / > 1 or > 1, 

Nevertheless, a specific orientation of anisotropy can be also chosen for some analysis. 



B. Case I: Propagation along the ^-axis of the PV 

PV-model-1: First, we consider a simple case where the excitation propagates only along the ^-axis in the PV. 
This type of propagation happens only when the propagation with the planar or convex wavefront reaches the facet of 
the PV in homogeneous media. For the sake of further simplicity, let the anisotropy be aligned in the circumferential 
direction such as d = d'^cf)^ > 1.0 in the block of anisotropy and = 1.0 everywhere. 

RA-analysis: The above conditions imply that A^i = 1.0, A^ = 0.0, d^ = 1.0 near the root of the PV, i.e.. 



at 6> 0. Then, by equation (29), we obtain the following equalities such as 

1 ^ If -cosO 



R2 \l-sm0^r/R 
and consequently the relative acceleration equation turns out to be 

d'^n' _ l-smO{l^r/R) fdO\ ^ -cosO dv 



dX^ (1 - sin (9 + r/R)^ \dn ^ 1 - sinO ^ r / R dn ' ^^^^ 

Note that each component depends on the angle 6 and the ratio r/R^ but its magnitude is relatively small and 
bounded. Moreover, the above conditions occur near the root of the PV, or ^ ^ implying 86 /dn ^ from n ^ (j), 
thus, by letting ^ = 0, the above equation boils down to 

9^n* _ R dv^ 

As expected, the above equation implies that there is no relative acceleration in the direction of the propagation^ but 
a certain magnitude of the relative acceleration can be generated in the direction of the wavefront caused by the 
geometry of the PV leading to the decrease of the propagational (conducting) velocity. 



Computational modelling: Figure 12 validates this analysis. For the planar propagation approaching the 
PV from the right wall, the propagational speed decreases in the PV due to the relative acceleration in the direction 
of the wavefront caused by the geometry of the PV. Since the propagation follows the direction of anisotropy in 
the atrium with increased conduction velocity, the decrease of the conduction speed can be more dramatic in the 
actual atrium as in vivo observed by Arora et. al. [4]. This phenomenon seems to confirm the dominance of the 
circumferential anisotropy in the actual PV as reported in ref. [44 . 

PV-model-2, RA-analysis and computational modelling: We also consider the ^-anisotropy aligned 
along the ^-axis such as d = d^O^ d^ > 1.0 in the block of anisotropy and d^ = 1.0 everywhere. But, for the sake of 
simplicity, we let d^ be constant along the (/)-axis. Since these conditions can be expressed such as 

dd^ ^ dd^ _ ^ dd^ _ dd^ dc/) dd^ dO _ dd^ dO 

dO dS dn 86 dn 86 dn 86 8n^ 
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with dO/dn « as before, we obtain the relative acceleration equation at the PV junction where ^ = such as 



d'^rf d{\ogdP)dv' 1 d'^d<^ 



R 



R 



dn 



(32) 



Observe that the first two components In the right hand side are the same as those in the relative acceleration 
equation (23) for anisotropy in the plane. The only additional factor is the last component which represents the 



relative acceleration caused by the unique shape and the anisotropy of the PV. In the last component, dilogdr) / dO 
has a nontrivial magnitude as shown in Figure [4|^, but is a negative value which is not summed up with the positive 
R/{R + r). Intuitively, this component indicates that the relative acceleration in the direction of the wavefront is 
compromised by the 6>-anisotropy that increases the conduction velocity. In results, the relative acceleration does not 



significantly increase. Figure [13] displays that the geometry or the 6>-anisotropy of the PV does not add significant 
amount of additional relative acceleration to yield conduction failure. 
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FIG. 12: PV with the 0- anisotropy 
(B), T = 850.0 (C). 



= 8.0 in (A). Initiated from the right wall, the membrane potential (u) at T — 750.0 
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FIG. 13: PV with the 6'-anisotropy 
(B), T = 750.0 (C). 



8.0 in (A). Initiated from the right wall, the membrane potential (u) at T = 650.0 
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FIG. 14: For 



4.0 and d = 7r/4, (aE/an)/E (A), G/E (B), and (aG/an)/E (C) versus the angle A. 
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FIG. 15: For ^ = (A), 6 — n/A (B), 6 — 11/2 (C), the sum of the components in equations (35) and (36) versus the angle A 



C. Case II: Propagation with an oblique angle to the PV 

PV-model-3: When the cardiac tissues are excited in the PV, the excitation starts to propagate along the ^-axis. 
Thus, we must take both of the (/)-direction and the ^-direction into account to understand the mechanism of the 
propagation in the PV. For this purpose, the relative acceleration approach can be used again, but the consideration of 
both directions without any simplification may bring incomprehensible complexity into the analysis, thus for the sake 
of simplicity, we suppose that the propagational direction and the direction of the wavefront are uniform along the 
wavefront. In addition, we also suppose that the diffusivity tensor dP and remain constant along the wavefront in 
the PV. These simplifications can be regarded to be reasonable considering the relatively small width of the myocardial 
tissue along the 6>-axis of the PV and the periodicity of and d^ along the ^-axis. However, they are just simpli- 
fications for mathematical analysis which may get its validation from computational modelling or clinical observations. 

RA-analysis: The above conditions lead to the following mathematical expressions for A^*, A^, and d such 
as 



dn 



dn 



0, 



dd^ 
dn 



dd^ 
dn 



0, for the wavefront in the PV. 



With the above equalities, the variable E and the derivative of E with respect to n are obtained as 

d^A? d^Al dE 2Rd^AlcosO dO 



E 



dE 
dn 



sin l9)+r)2' dn {R(l - smO) ^ rf dn' 
Similarly, the variable G and the derivative of G with respect to the parameter n are obtained such as, 

d^AecosO dG _ d^Ae{smO{R^r)-R) f dO^ 



G 



dG 
dn 



-sin l9) +r' dn - sin l9) + r)^ \dn 

If we let the propagational direction be at the angle A with respect to the ^-axis such as 

Aq — cos^, A^ — sin^, for < ^ < tt. 



(33) 



(34) 



Then, we get dO/dn = — sin^ and equations (33) and (34) are rewritten as 

d'f' sin^ A 



E = 



d^ cos^ A 



dE 
dn 



2i?d^ cos(9sinM 



and 



R^ -sin l9) + r)2' dn - sin l9) + r)^ ' 

d' cos cos A dG d^{{R-\-r)smO — R)cosAsmA 



-sin i9) +r' dn 



- sin l9) +r)2 



(35) 



(36) 



Now, each component is the function of {d^^d^) for anisotropy, (i?, r) for the geometry of the PV, and {O^A) for 
the angle between anisotropy in the geometry and the propagational direction. After being divided by E, the above 
components are plotted in Figure 14 where we set d^ = d'^ = 4.0 and = 7r/4. We observe that each component has 



different magnitude and sign depending on the angle A. Consequently, the summation of the components also displays 
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the strong dependence on the angle A and for diverse relative acceleration as shown in Figure 15 In words, this 
means that the relative acceleration in the PV with anisotropy is rapidly changing depending on the angle between 
the propagational direction and anisotropy or the geometry of the PV. Note that the sensitivity to the angle A is 
significantly large and rapidly varying compared to that of anisotropy in the plane as shown in equation (24). 

For example, as shown in Figure 15 A. for ^ = and Figure fl5b for ^ = 7r/4, the relative acceleration is large when 



the propagational direction is at an angle of 7r/4 with the ^-axis, but the relative acceleration is negligible when it 
is close to 37r/4. Interestingly, ^ = and A = tt yield the similar amount of the relative acceleration, but in the 
different signs. Accordingly, whether the propagation direction is toward the plane (i.e., toward the atrium in real 
anatomy) or away from the plane may play a critical role for conduction failure. However, regardless of the angle A, 
the relative acceleration becomes weaker as is closer to 7r/2, i.e., far from the atrium, thus no significant conduction 
failure is likely to occur (Figure fTsfc) in the area of the PV with ~ 7r/2. 



Computational modelling: To demonstrate the validity of the above analysis, a T-shaped domain is de- 
signed by placing a PV in the middle of it. The left rectangle is — 150 < < and —70 < z < 70, while the right 
rectangle is —0 < y < 150 and —150 < z < 150. At the center of (0,0,0), a PV is constructed with R = 70 and 
r = 35. Moreover, anisotropy is aligned with = dP = 8.0 in the area of —30 < ^ < of the PV which is defined 
exclusively in the range of < x < 60. 

With a point-initialization in the right rectangular domain, the only pathway to propagation to the left rectan- 
gular domain is through the PV with the ^-anisotropy and the (/)- anisotropy. When the excitation is initiated at 
(0.0, 100.0, 0.0) as shown in Figure 16 A., the anisotropy of the PV does not generate a sufficient relative acceleration 
to block the propagation, since the angle A seems to be close to 7r/2. As a consequence, the excitation propagates 
through the P V without conduction failure . See Figure [l6b and p!6t 



On the other hand, when the excitation is initiated at (0.0, 100.0, 80.0) as shown in Figure 17 A., the excitation 



propagates toward the PV with an oblique angle A which is more acute than the previous case and is possible around 
7r/4. As displayed in Figure 15, this acute angle A induces a large relative acceleration to cause conduction failure as 
shown in Figure [iTp. Since the only difference between the two simulations is the location of point-initialization and 
subsequently the angle A toward the anisotropy in the PV, these models seem to confirm the previous analysis from 
the relative acceleration equation for the PV with anisotropy. 
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FIG. 16: PV with anisotropy of = < 
potential {u) at T = 800.0 (B) and T 



^ = 8.0 and the location of point initialization at (0.0, 
1000.0 (C). 



100.0, 0.0) (A). The membrane 




FIG. 17: PV with anisotropy of (f ^ ^ 8.0 and the location of point initialization at (0.0, 100.0, 80.0) (A). The membrane 
potential {u) at T = 600.0 (B), T = 850.0 (C). 
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IX. DISCUSSION 



The simplicity of the relative acceleration approach is gained from the use of a tensor which measures the 
spreading-out or convergence of the trajectories to indicate the changes of the propagational speed along the wavefront. 
The unification of the geometry and anisotropy also contribute to this simplicity of analysis in view of the Riemannian 
geometry. The six problems on simple surfaces and three models of the PV are analysed by the relative acceleration 
approach and compared to the computational simulations and clinical observations which seems to provide the validity 
of the proposed approach based on the proposed hypothesis. For a practical application, the last model of the PV 
which shows the strong angular dependency of the PV for cardiac excitation propagation may provide a model for the 
mechanism of unidirectional pathway [33 , a directional-sensitive pathway of cardiac excitation propagation, which is 
regarded as a necessity for atrial reentry resulting in atrial fibrillation. 

In spite of these successful examples, there are currently some restrictions of applying this approach directly for 
real cardiac electrophysiological phenomena. The first restriction is due to the assumption on the traveling wave 
solution for cardiac excitation propagation in the homogeneity of the excitable media. Therefore, the shape of the 
cardiac action potential should be roughly maintained in the propagation for valid analysis. For example, the dramatic 
upstrokes or downstrokes, the rapid slow-down or speed-up caused by the heterogeneity of the cardiac tissues cannot 
be dealt with the relative acceleration approach unless other approximations are newly introduced. In fact, anisotropy 
also induces the changes of the shapes of the cardiac action potential, but the amount of it is relatively small, thus the 
analysis from the relative acceleration approach is not seriously affected. The second restriction is regarding on the 
singularities of anisotropy. The relative acceleration equation is derived in the neighbourhood of a point in a curved 
surface, thus the discontinuity of anisotropy is not a problem, but the singularity on the point is. For example, when 
multiple anisotropics come across a point, it is equivalent to a singularity point in mathematics and this case is not 
considered in this paper. 

Nevertheless, the relative acceleration approach is an attractive option because it can be easily extended to more 
complex higher-dimensional anisotropic surfaces such as multiple layers of curved surfaces or three dimensional 
anisotropic space which are prevalent in the real heart. Also, revealing the relationship between the kinematics 
approach and the relative acceleration approach, or in other words the relationship between the shape of the wave- 
front and the distribution of the trajectories, may yield more practical tools applicable to clinical and surgical planning. 
In view of physiology, one of the most important consequences of this approach is the use of the trajectories and its 
validity in the phenomena of cardiac excitation propagation. This result seems to implicitly suggest the existence of 
the moving bodies^ or at least something equivalent, in cardiac excitation propagation and it seems to open questions 
on what is actually moving along the trajectories in cardiac excitation propagation. 
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In the following calculations, the sum is taken over all the indices repeated above and below except k. Let's 
differentiate ^/gg^^ with respect to the axis Xk- Using the chain rule and using d-y/g/dxk = ^/g^^k^ obtain 
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X. APPENDIX 



A. Differentiation of y^g' 



kk 




(37) 
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To replace the differentiation of g^^ with the differentiation of the metric tensor gj^j^, we use the identity of g'^^g^j = 
to obtain 



and the equation obtained by differentiating g^^gt^i — as 



(38) 



dxk oxk 



Note that this reduce equation ( |38[ ) to the fohowing: 

v^5'''r;:,-V^/V'^. (39) 



Replacing the differentiation of g^^ with the equalities in r^j/^ and F*^, which are Christoffel symbols of the first kind 
and the second kind respectively, such that 



With these equalities and the condition of orthogonality of the curved axis x^^ equation (39) becomes 

1 5(V^5"=^ 



^/9 dXk 



5"(0-ry, m^k (40) 



For surfaces of revolution such as spherical shell and torus, the curved axis Xm representing the rotational direction 
is in the direction of the Killing vector, or the Killing direction^ which satisfies the condition dg^jy/dxm = j30] . 
For example, the Killing direction is the axis for spherical shell, the azimuthal angle, and the (f) axis for torus, the 
rotational direction of the radius of a great circle. If the axis x^ is in the Killing direction, then the above equation 
is non trivial only when 

I d ^ m^fc.D (41) 



^/g dXk 



B. Geometric factors of the PV Column 



R : radius of the circular arc, r : radius of the circular hole, 
6 : longitudinal angle from the root of the column, (j) : circumferential angle , 



Parameterization and its differentiation 



x= {R{1 - cos l9), - sin l9) + r) sin0, {R{1 - sin l9) + r) cos^), 

X6> = {R sin 9^—R cos sin (/), —R cos 9 cos (j)) , 

X0 = (0, {R{1 - sin6>) + r) cos(/), -{R{1 - sin6>) + r) sine/)). 



Metric tensors 



gee = • = R^, g^^ = x^ • x^ = {R{1 - sin6>) + r)^, y/g = R{R{1 - sinO) + r), 

H ^9ee_ ^ 1 ee ^ ^ 

^ g (i?(l-sin^) + r)2' ^ g R^' 
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Christoffel symbols and its differentiation 



1 



RcosO{R{l-smO)-\-r), 



1 



2 dO 



= Q^'^BH = cos^((l - sin^) + r/R), T^^ = g^^T^^e 



-RcosO{R{l-smO)-\-r), 
— cos ^ 



(1 - sin (9) +r/i?' 



do 



{R{1 - sin (9) + r)2 ' 00 {R{1 - sin (9) + r)3 ' 
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